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Abstract 


A computationally efficient 3DHZETRN code capable of simulating High Charge (Z) and Energy 
(HZE) and light ions (including neutrons) under space-like boundary conditions with enhanced 
neutron and light ion propagation was recently developed for a simple homogeneous shield object. 

Monte Carlo benchmarks were used to verify the methodology in slab and spherical geometry, and the 
3D corrections were shown to provide significant improvement over the straight-ahead approximation 
in some cases. In the present report, the new algorithms with well-defined convergence criteria are 
extended to inhomogeneous media within a shielded tissue slab and a shielded tissue sphere and tested 
against Monte Carlo simulation to verify the solution methods. The 3D corrections are again found to 
more accurately describe the neutron and light ion fiuence spectra as compared to the straight-ahead 
approximation. These computationally efficient methods provide a basis for software capable of space 
shield analysis and optimization. 

Introduction 

Early space radiation shield code development relied on Monte Carlo (MC) methods [Alsmiller 1967, 
Lambiotte et al. 1971] and made important contributions to the space program. Due to intensive computational 
requirements, MC methods utilized restricted one dimensional problems leading to imperfect representation of 
appropriate boundary conditions [Alsmiller et al. 1972, Pinsky et al. 2001, Armstrong and Colburn 2001, Foelsche et 
al. 1974]. Even so, intensive computational requirements for MC codes remained, and shield evaluation was made 
near the end of the design process in greatly simplified geometry to enhance computer efficiency [Armstrong and 
Colburn 2001, Wilson et al. 2002]. Resolving shielding issues at the end of the design cycle had a negative impact 
on the design, since resolving issues early could have minimized shield augmentation requirements and associated 
launch costs. This is especially true in post-launch augmentation, as was done for the International Space Station 
(ISS) [Shavers et al. 2004]. Furthermore, added shielding at the end of the design process could require de-scoping 
mission objectives, as instruments and equipment may be removed to meet launch requirements, as was done on the 
Viking Project. 

Improved spacecraft shield design requires early entry of radiation constraints into the design process to 
maximize performance and minimize costs. As a result, NASA has been investigating high-speed computational 
procedures to allow shield analysis to be part of the preliminary design concepts following through to the final 
design allowing shield optimization procedures [Wilson et al. 2003a, 2004a,b]. For the last several decades, NASA 
has pursued deterministic solutions of the Boltzmann equation allowing field mapping within the ISS in tens of 
minutes [Wilson et al. 2007] using standard finite element method geometry common to modern engineering design 
practice [Qualls et al. 2001, Wilson et al. 2003a, 2004a]. Ray tracing procedures in complicated geometry hinders 
the application of MC methods to such engineering models. Even so, deterministic methods have relied on the 
straight-ahead approximation, resulting in the HZETRN code with loosely defined impact on model uncertainty 
[Wilson and Khandelwal 1974] but it handily facilitated the ISS augmented design for which no other code was 
capable. Recently, a 3D version of HZETRN has been developed in simple geometry (homogeneous sphere) with 
improved convergence criteria and verification of 3D effects using MC methods [Wilson et al. 2014]. 

Herein, the homogeneous media limitation will be removed from the solution methodology, and MC codes 
are used again to verify the solution methods in simple geometry. The restriction to slab and spherical objects will 
be maintained to allow efficient MC simulations for verification purposes. The next step in this development will be 
to use more complex geometric models so that simple mapping of the present methodology into more realistic 
applications can be studied. In the present report, the current status of transport code development will be briefly 
reviewed with emphasis on extending these developments into a generalized 3D version of the HZETRN code. 
These advances will use available MC codes, Geant4 [Agostinelli et al. 2003], FLUKA [Fasso et al. 2005, Battistoni 
et al. 2007], and PHITS [Sato et al. 2006, 2013] to judge the veracity of these developments, especially with regard 
to their 3D aspects. 

Deterministic Code Development 

The relevant transport equations are the linear Boltzmann equations derived on the basis of conservation 
principles [Wilson et al. 1991] for the flux (or fiuence) density, , of a j type particle in the continuous 
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slowing down approximation (CSDA) in which atomic processes are described by the stopping power Sj(E) for each 
ion type j (vanishes for neutrons, 7 = n) as 



where the differential operator on the left hand side is defined as 

B {x, n,E)] = n. Vcf>^ (x, {E)4>^ (x, i?) ] + {E)cf>^ (x, El, E) , (2) 

and solved subject to a boundary condition over the enclosure of the solution domain as shown in Fig. 1. 



Fig. 1. Geometric relations of quantities useful in solving equation (1). The symbol n is a unit normal vector. 


Application of the CSDA in both laboratory and space shielding has been wide and the resulting errors 
discussed elsewhere [Wilson et al. 1984, 1987a, 1994, Shinn et al. 1998, Tweed et al. 2006a, 2006b]. Equation (1) 
can be rewritten as a Fredholm equation [Wilson 1977] given by 
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In equations (1) and (3), F(Q,x) is the vector connected to x along -Q; p is the projection ofx onto (Fig. 1), and d 
is the projection of F(Q,x) onto fl. The quantity Rj(IE) is the distance a type j ion with energy E will travel before 
losing all of its energy to excitation/ionization of atomic electrons, and Pj{IP} is the probability a type j ion with 
energy E will not have a nuclear interaction in coming to rest in the media. The zero charge limit of equations (1) 
and (3) is discussed elsewhere [Wilson et al. 1991]. The range-energy relation is given by 
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and the nuclear attenuation function is given by 
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where Oj{E ') includes nuclear elastic and reactive processes. One obstacle to solving either equation (1) or (3) is the 
need to evaluate the integral dCV at arbitrary locations within the media and development of computational methods 
to efficiently handle this limitation. 

Atomic interactions limit the contributions of charged particles in the transport process. For example, the 
protons and alphas produced in aluminum below 100 MeV/n contribute to the fluence only within a few centimeters 
of their collisional source, and the heavier ions are even more restricted. This is an important factor in that the 
transported secondary charged particle flux tends to be small at low energies, and the role of additional nuclear 
reactions by these low energy ions are likewise limited, an important fact that is implemented in the current 
development. 

As noted above, a prime limitation in evaluating equation (1) or (3) is the evaluation of the integral over 
do.' at arbitrary locations within the media. The approach to a practical solution of equation (1) or (3) is to develop a 
progression of solutions from the simple to increasingly complex allowing early implementation of high- 
performance computational procedures and establishing a converging sequence of approximations with established 
accuracy criteria and means of verification and validation. The first step leading to the lowest order solution reduces 
the evaluation by introducing the straight-ahead approximation as guided by the nucleon transport studies of 
Alsmiller et al. [1965] using MC methods in which the differential cross sections were approximated as 


a^YE.E'MM') = <J.,{E,E')6{^1 - O') . (8) 

This approximation produced dose and dose equivalent results to be within the statistical uncertainty of the 
MC result obtained using the fully angle dependent cross sections in slab geometry [Alsmiller et al. 1965, Wilson et 
al. 1991]. The relation of angular dependent cross sections to spacecraft geometry was further examined by Wilson 
and Kandelwal [1974] using asymptotic expansions about angular divergence parameters, demonstrating errors in 
the straight-ahead approximation to be on the order of the square of the ratio of distance of divergence (few 
centimeters) to radius of curvature of the shield (few to several m), resulting in a small relative error ~10 "' in most 
human rated space systems [Wilson and Khandelwal 1974]. Introducing the straight-ahead approximation, equation 
(8), into equation (3) reduces the Fredholm equation to a Volterra equation that can be solved using marching 
procedures [Lamkin 1974, Wilson and Lamkin 1974, 1975, Wilson 1977, Wilson and Badavi 1986]. Considering the 
success of the straight-ahead approximation for nucleons encourages its use in HZE transport where the greater mass 
of heavier ions reduces even further the angle of scatter [Wilson 1977]. The success of equation (8) results from the 
fact that lateral diffusion from a given ray is compensated in a flat plate by lateral diffusion along adjacent rays. 

Numerical marching procedures were developed to solve the transport equation under the straight-ahead 
approximation, resulting in the HZETRN code. A corresponding nuclear fragmentation model, NUCFRG2, was also 
developed for HZETRN, and the verification and validation processes utilizing the NUCFRG2 database are 
described elsewhere [Wilson et al. 1987b, 2005, 2006]. 

Space flight validation of HZETRN has been limited almost exclusively to low Earth orbit (LEO), 
containing both trapped particle and attenuated galactic cosmic ray (GCR) components. The two primary limitations 
in the LEO trapped environmental models AP8MIN and AP8MAX as discussed by Wilson et al. [2003b] is the 
assumption that the trapped particles are isotropic (resulting from the omnidirectional fluence description) and poor 
representation of the dynamic behavior. These omnidirectional models have been relatively successful in describing 
the radiation environment aboard the highly maneuverable shuttle spacecraft wherein anisotropies tend to be 
averaged. Such models have been found to be less accurate in the formation flying of the ISS, mainly oriented in the 
local horizontal plane along the velocity vector as was demonstrated elsewhere [Hugger et al. 2003; Wilson et al. 
2005; Nealy et al. 2006; Slaba et al. 2011a, 2013]. A dynamic and anisotropic trapped proton environmental model 
was developed for future use in LEO shield design and operations [Wilson et al. 2003b]. More recently, an updated 


3 


trapped proton model, AP9, has been made available and preliminary validation eomparisons have shown reduced 
uncertainties [Badavi et al. 2014]. 

These environmental models are placed in a suitable form for evaluation of the incident radiation on the 
bounding surface of the six-degree of freedom motion of an orbiting spacecraft for shield evaluation [Wilson et al. 
2005]. To test the dynamic behavior, shuttle thermoluminescent dosimeter (TLD) data from 1983 to 2000 were 
used, giving good coverage for nearly two solar cycles. A sample of shuttle TLD measurements and the 
computational model results is given in Table 1. With the use of a normalization procedure intended to represent 
environmental model uncertainty, Badhwar et al. [1996, 2001] showed that the root mean square error of both 
observed and calculated dose rate and dose equivalent rate is within 15%. The normalization procedure is based on 
scaling the model trapped proton contribution to the TLD measurements. Recent validation studies by Slaba et al. 
[2011a, 2013] did not utilize normalization procedures and focused entirely on the GCR component of the LEO 
environment. Those studies revealed larger differences exceeding 40% if pion production and the associated 
electromagnetic cascade were not included, as was the case in most of the previous validation efforts. These 
comparisons highlight the usefulness of HZETRN within the straight-ahead approximation in many cases, but that 
further improvement is still necessary. 


Table 1. Comparison of straight-ahead HZETRN transport model with shuttle flight data. [Wilson et al. 2003b]. 


Flight 

Date 

DRNM* 

DLOC 

TLDt (pGy/d) 

Calc. (|iGy/d) 

STS-41A 

11/83 

6421 

3 

64.6 

59.6 

STS-51D 

4/85 

6661 

4 

917.4 

889.3 

STS-31 

4/90 

5701 

1 

2141 

2290 

STS-43 

8/91 

5894 

4 

20.7 

18.6 

STS-62 

3/94 

6771 

1 

94.3 

89.2 

STS-65 

7/94 

6822 

2 

28.3 

25.1 

STS-67 

3/95 

6925 

3 

250.8 

238.1 

STS-80 

11/96 

6973 

4 

264.4 

256.5 

STS-82 

2/97 

7074 

1 

2978 

3080 

STS-91 

6/98 

6894 

1 

89.1 

83.2 

STS-101 

5/00 

6460 

2 

140.8 

131.1 

STS-92 

10/00 

6417 

2 

165.9 

153.4 


* Deep River Neutron Monitor count rate, f GCR corrected TLD- 1 00 data 


Most of the previous validation calculations were performed in a homogeneous media (aluminum alloy 
2219), but radiation protection requires, at a minimum, the calculation of the environment within the human body, 
yielding an inherently inhomogeneous problem. 

Applying radiation exposure constraints to the space radiation shield design process requires not only the 
estimation of the internal spacecraft environment but how that environment is transported in human tissues of a 
characteristic human shape having location and orientation within the spacecraft [Wilson et al. 1993]. To examine 
this problem experimentally, a human skeleton (head and torso) was covered and fdled with a tissue equivalent 
plastic and sectioned to allow instrumentation for in-flight measurements (see Fig. 2). The torso phantom was 
Computerized Tomography (CT) scanned to define the geometry for use as a computational model, especially in 
regard to organ location. The instrumentation used TLD, nuclear track foils and Tissue Equivalent Proportional 
Counters (TEPCs) to allow the mapping of the dose and dose equivalent within the phantom torso model. The 
instrumented torso phantom was flown on STS-91 followed by data analysis and dosimetric analysis using the 
HZETRN transport code. The results of the measurements expressed herein as the average quality factor are given in 
Table 2. If the environmental model is again normalized to measurements, the calculated dose equivalent was within 
12% of the measurements except for the stomach which was 19%. The cause of the differences was in large part 
due to the need to correct the galactic cosmic ray under-response of the TLDs [Cucinotta 2004], and the corrected 
data were all within 12% of the computational models. 

The usefulness of the TEPC in evaluating dose equivalent in the space environment being established, 
Badhwar provided a test rig of five TEPCs (Fig. 3) to study dose equivalent attenuation in the LEO environment 
within the shuttle spacecraft. The bare detector allows characterization of the ambient field within the shuttle while 
the shielded detectors (spheres in the figure) sample the modified environment due to the shielding material. A 
study of the LEO radiation attenuation in polyethylene was conducted in the shuttle bay on STS-81 and are 
compared with a similar study of aluminum shielding on STS-89. 
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Spacecraft power 

Power distributor 
and switch box 


Bare | 
detector 



8 in. 

sphere sphere 


Fig. 3. Apparatus for shield attenuation studies during shuttle flight [Cucinotta 2004]. 


Table 2. Calculated and “measured” (STS-91) quality factors. All dose equivalents were within 12% except for 
stomach which was -19.7%. 


Organ 

Qcq\c 

^meas 

Brain 

2.10 

1.8 0.1 

Colon 

2.17 

2.0 0.2 

Heart 

2.15 

2.0 0.4 

Stomach 

2.20 

1.7 0.3 

Thyroid 

1.89 

1.7 ±0.2 

Skillbreast 

2.10 

1.8 ±0.2 

Skillabdomen 

2.11 

1.9 ±0.2 


Note that a direct comparison of aluminum shielding effectiveness to that of polyethylene cannot be made 
as the two flights are in radically different environments. STS-81 is in a high inclination orbit (51.6°) with large 
contributions from galactic cosmic radiation and STS-89 is in a low inclination orbit (28.5°) where the environment 
is dominated by trapped protons. The results are compared with the HZETRN model calculations in Fig. 4. 

One TEPC failed on STS-89, but it is still clear that good agreement was found between the FIZETRN and 
the measured attenuation curves. This is also true for the STS-81 results for polyethylene shields as shown in the 
figure. It is clear that polyethylene has superior attenuation characteristics even in environments with large 
contributions from galactic cosmic radiation (STS-81). As a result, polyethylene is the material of choice for 
augmenting the ISS shield in spite of its need to be encapsulated to guard against fire hazards. Most important is the 
ability of HZETRN, even in the straight-ahead approximation, to evaluate correctly the dose equivalent from all 
components important to human rated systems and provide estimates of ISS augmentation requirements and the 
value of the benefits. Consequently, the layout of the augmentation in the ISS was evaluated using the HZETRN 
code within the straight-ahead approximation [Shavers et al. 2004]. 
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Fig. 4. Experimentally measured and calculated attenuation curves in aluminum and polyethylene in shuttle flight 

[Cucinotta 2004]. 


It is clear that the straight-ahead approximation provides accurate dosimetric results and a good 
approximation to the particle fluence in many circumstances (see Badhwar et al. [1995]). One limitation of the 
straight-ahead approximation is near the front boundary where the calculated neutron fluence vanishes (unless 
neutrons are part of the external radiation environment). This results from the straight-ahead approximation and can 
be improved by employing 3D corrections containing backward leakage to the front boundary. While leakage is an 
important process of neutron transport in space systems, it is less important to the charged particle transport due to 
range/energy relations. 

In the early transport studies of Wilson and Lamkin [1974, 1975] and Lamkin [1974], it was demonstrated 
that once neutrons are produced by a proton beam, that the recoupling of neutrons back to the proton field is limited, 
but the neutron fields tend to build-up with increasing penetration depth. This mainly occurs since many protons 
produced by neutron-induced reactions are of lower energy than the forward components and are quickly removed 
by atomic interactions (especially the isotropically produced protons). The probability of even a 100 MeV proton 
undergoing a nuclear reaction is only several percent, and the probability decreases substantially at lower energies. 
The low-energy neutrons, on the other hand, have no effective atomic interactions and propagate until a nuclear 
interaction occurs or they escape from the media. Yet, these low-energy protons and other light ions produced in 
tissue through neutron collisions contribute significantly to biological injury and must be accounted. This was the 
basis of the studies of Heinbockel et al. [2003] starting with Clowdsley et al. [2000, 2001] using the NUCFRG2 
database [Wilson et al. 1995, 1998], with recent improvements provided by Slaba et al. [2010a]. Herein, we will 
concentrate on extending the transport formalism to full 3D using the NUCFRG3 [Adamczyk et al. 2012] database. 
Future work will be a direct evaluation of equation (3) quantifying the propagated error of the 3D neutron code. 

Extension of the transport formalism begins by writing the nuclear reactive and scattering differential cross 
sections in the following form 




(9) 


where the first term is isotropic and associated with lower-energy particles produced, including target fragments, 
defined by 


= 2 J^^^a^,{E,E\n,n')dn\ (10) 

where 2nB represents the backward hemisphere, and the second term in equation (9) is highly peaked in the forward 
direction and is associated mainly with direct quasi-elastic events and projectile fragmentation products [Wilson 
1977, Wilson et al. 1988]. In this context, quasi-elastic refers to collisions in which the energy transfer is small 
compared to the incident energy, and the secondary products are highly peaked in the forward direction. 
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Surprisingly, even nucleon-induced reactions follow the simple form of equation (9), and the isotropic term 
extends to relatively high energies (see Fig. 5). Further details on the nature of the elastic and reactive cross sections 
used in HZETRN are given in Wilson et al [1991] and Cucinotta et al. [1994, 1998]. 



Fig. 5. Isotropic (iso) and forward (for) neutron spectra produced by 500 MeV protons in aluminum. 


The first step towards a more general solution uses a bi-directional neutron transport approximation as first 
investigated by Clowdsley et al. [2000, 2001] using multigroup methods and recoupling to the light ions using an 
analytic solution as a quadrature over the induced light ion source term. These types of approximations were further 
studied and improved by Slaba et al. [2010a] and Heinbockel et al. [2011a, 2011b]. 

In the most recent bi-directional model described by Slaba et al. [2010a], the charged particle components 
were first solved in the straight-ahead approximation with the neutron coupling represented only by the forward 
component. The isotropically produced neutrons are then solved using a coupled bi-directional transport formalism 
to complete the neutron solution [Slaba et al. 2010a]. These isotropically produced neutrons further collide with the 
media, giving rise to an isotropic source of ions that are transported assuming only atomic interactions as an 
approximation to the complete solution, justified by the limited range of these ions. The governing equation is given 
by equation (1) which is solved using equation (8) for charged components and assuming the separation of the 
neutron field into a forward and isotropic component as in equation (9). The charged components are solved using 
the straight-ahead approximation as 


E/; ' ^,kiE.E')4^kjorME')dE' , (11) 

for j ^n, and for j = n, only the forward produced neutrons, given by equation (9), are represented as 

[n.V + aSE)]<j>^j,AxM.E) = jy^,j„{E,E')c^ujoAxM,E')dE'. (12) 

Equations (1 1) and (12) are solved simultaneously using methods previously developed [Wilson et al. 2006; Slaba et 
al. 2010b]. Note that the ions generated by the isotropically produced neutrons are neglected in the above but will be 
treated later. Whereas the full straight-ahead cross section of equation (8) appears in the ion transport of equation 
(11), only the forward components of the cross sections appear in the transport of neutrons in equation (12) given by 

<^nkjor{E.E') = (13) 


The isotropic source of neutrons are treated as a perturbation given by 
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(14) 


. V + (x, Q, E) = /“ {xM\E')<m'dE' 

+E XT X. V {xM\E')dn'dE'. 

k 

Couplings between the isotropic neutron field and the ion fields of equation (11) are not explicitly included in 
equation (14) but are treated once the isotropically produced neutrons are evaluated. Clearly, equations (11) and (12) 
can be solved independent of equation (14), but the forward solution of equations (11) and (12) enters into equation 
(14) generating the isotropic neutron source. Whereas equations (11) and (12) take on value only if an external 
source is present at the boundary, equation (14) has no external sources, and the resulting field develops from the 
internal isotropic neutron source. Development of an algorithm for solving equation (14) is expedited by 
approximating the double differential cross sections as 

a^,^{E,E\n,n') = cr,^^^j{E,E')6{n - n') + a^^,{E,E')6{n + n') , (15) 


where 


<^nuAE,E') = f^^^a^JE,E',f2,f}’)dii' = a^^j„{E,E') + 2na^„^^JE,E') / dvr. 


(16) 


with 2nF as integration over the forward hemisphere and 


^ nn.h V 


{E.E') = f^^^a^^(E,E\n,n')dn' = 27ra„,^^^JE,E') / 4^, 


(17) 


with InB as integration over the backward hemisphere. Similar relations for a^f, in equation (14) are given by 

^uk,.JE,E\n,n') = - n') + ( 18 ) 

where 

<^nk,Mf)^E,E') = J^^^a^,aE,E\n,n')dn' = 2na^,^^jE,E') / 4^, (19) 

with 2nF as integration over the forward hemisphere and 

<^nk,so(,){E,E^) = = 2na^,^^JE,E') / dvr, (20) 

with 2nB as integration over the backward hemisphere. Using equations (15) to (20) in equation (14) results in 

[n.V + a^{E)]^^^^Jx,n,E) = J"a^^j{E,E')<l^^^,Jx,n,E')dE' 

+J^^nu,>,(^^E')4>^^^Jx,-n,E')dE’ 

+E XX ^ ME')dE' 

k 

+E XX ^ '^hfor (x,-^,E')dE' 
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(21) 


where the backward moving isotropic induced neutron field 4>^ satisfies a similar equation obtained by 

replacing Q by -Q in equation (21) as 


J »oo 

^ <^nnAE,E')cj^^,soi.^-n,E')dE' 

J »oo 

, ^un,,{E,E')4>^^^Jx,n,E')dE' 

+E /“ ^n,,soU) (E, E ')<!>, {X, -n,E')dE' 

k 

+E J^^n,,soii)iE,E')cj^,jJx,^l,E')dE' 


( 22 ) 


Note that (j)k for ™ equations (21) and (22) represents the forward directed flux along ft, and consequently, 

for (*’ E') = 0 . The isotropically produced neutrons induce a source of isotropic light ions in collisions with 

the media resulting in a light ion fluence given by 


n.V(^^{xM,E) - ^j-[s^{E)c^^{xM,E)\ cl,,^,Jx,n,E) = a^jE,E')cl,„^^Jx,n,E')dE\ (23) 


where j #«, and (a;, JT, i? ') is the isotropically produced neutron fluence. Equations (21) and (22) are solved 

simultaneously for no external fluence incident on the media directed along Q (see Fig. 1). Note that the implicit 
assumption of equations (21) and (22) is that lateral diffusion from the ray is fully compensated by diffusion from 
adjacent rays. This is strictly true for a semi-infinite slab and represents the use of periodic boundary conditions as 
an importance sampling method common to MC simulation practice. Note that the structure of these equations is the 
same as described by Wilson et al. [2014] for uniform media in simple geometry, and herein, the continuity of 
solution across the interface of the two media is enforced in the numerical solution. The formalism described by 
equations (1 1) - (23) is hereafter referred to interchangeably as the bi-directional solution, or the A^= 2 solution. The 
N ~ 2 terminology was discussed by Wilson et al. [2014] and indicates that the isotropic neutron field is described 
by two ray directions, forward and backward. 

The solution is now examined near the interfaces of a 30 g/cm^ semi-infinite slab of ICRU tissue (mass 
composition of 76.2% O, 11.1% C, 10.1% H, 2.6% N) [ICRU 1993] shielded and backed by 20 g/cm^ of aluminum 
exposed to the Webber solar particle event (SPE) model [Webber 1966] given as 

^(.E) = 12!fc^exp{[239.1 - p{E)] / lOO} , (24) 


with p{E) = ^ E{E + 1876) . Transport procedures were carried out using the bi-directional {N = 2) solution with 
results shown in Fig. 6. 

Results are shown for the usual biologically significant depths assigned by the ICRU [1993] as 0.007, 0.3, 
and 1 g/cm^ but include fluence in the aluminum shield adjacent to the forward (near) and the backward (distal) 
interface. One first observes that the neutron fluence varies little across both regions. The proton and alpha fluences 
are nearly constant in the aluminum shield before the near interface but vary rapidly within the tissue after passing 
the interface and approach a new equilibrium spectrum by 1 g/cm^. The main change is in greater thermalization in 
tissue. The proton and alpha fluences in tissue before the distal interface are nearly constant but vary rapidly in the 
aluminum past the interface. 

The solutions for charged components at low energies are proportional to the inverse stopping power of the 
media, and the transition of the charged particle spectra in the first several hundred microns of the interface is clear 
from the figure. Note that no transition occurs in the aluminum before the near interface and in tissue before the 
distal interface resulting from the use of the straight-ahead approximation for charged components. 
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Fig. 6. Particle fluence induced by the Webber SPE in 30 g/cm^ ICRU tissue slab shielded and backed by 20 g/cm^ 
aluminum slabs in the vicinity of the forward (left pane) and backward (right pane) interface. The depths of 0 and 30 
correspond to the front and back interfaces, respectively. Fluences are shown at the biologically significant depths 
(solid) in tissue defined by ICRU 51 [1993] and similar depths (dashed) in the aluminum shield near the interfaces. 

An expanded view of the alpha fluence on a linear scale is given in Fig. 7 for several depths in the 
neighborhood of both the near and distal interfaces. It is clear that a large transition occurs as a result of the change 
in both stopping power and alpha production within the two media. Improvements in the transition description will 
require, at a minimum, a bi-directional treatment of the charged particle fluence. 

This bi-directional approximation {N = 2) in slab geometry is compared with the MC codes, Geant4, 
PHITS, and FLUKA, in Fig. 8. All Monte Carlo fluence results presented herein are plotted as mean values with 
estimated statistical uncertainty represented by one standard deviation and shown as vertical error bars placed at the 
energy bin midpoint. Details of the MC simulation setups are given in Appendix A for Figs. 8 and 17-23. Although 
the bi-directional approximation makes significant improvements to the neutron fluence in slab geometry, additional 
improvements with further 3D corrections are expected, especially for the charged particle fluences near interfaces. 
The success of the straight-ahead approximation under validation in human rated systems (shuttle, ISS) clearly 
meets the requirements for a useful solution of the transport problem, although the bi-directional solution of the 
neutron fluence makes marked improvements. Considered in the subsequent section is the transport in finite 
inhomogeneous objects in which lateral leakage can play an important role. 




Fig. 7. Alpha fluences from Fig. 6 with the vertical axis on a linear scale and the horizontal axis restricted below 10 

MeV/n. 
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Fig. 8. Comparison of fluence spectra at 0 g/cm^ (left pane) and 30 g/cm^ (right pane) in a 30 g/cm^ ICRU tissue slab 
shielded and backed by 20 g/cm^ aluminum slabs exposed to the Webber SPE spectrum. 


3D Marching Procedures 

A sequence of approximations that span the formalism from the simple straight-ahead approximation to 
more complex propagation algorithms is described in this section. The governing equations in the CSDA ignoring 
Coulomb scattering are given by equations (1) and (3), and the cross sections include nuclear elastic and reaction 
processes. The double differential cross sections are found to be closely approximated by a forward component and 
an isotropic component given in equation (9), where the isotropic cross section is given by equation (10). The 
corresponding forward component is found from 


As was done in the development of the bi-directional transport formalism, the solution is divided into forward and 
isotropic components where the isotropic neutron term is treated as a perturbation given by the following equations 
(26) - (29) 




(26) 


On the scale of a vehicle or habitat, the space environment is nearly uniform in position and varies slowly in angle 
(if not isotropic) and each directional component Qo (see Fig. 9) is propagated according to equation (26) or 
according to the corresponding Volterra equation as 


S.,(E)P,(E) 


■ r(r^j,a:),S7jpi? 


A. 




(27) 


+ - 


SAE)R{E)^Je Je 


where = x + [Rj(E) — R-{E')]Ttg. The means of handling the zero neutron charge is discussed elsewhere 
[Wilson etal. 1991]. 
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Fig. 9. Geometry of 3D marching procedure. 


Whereas the full straight-ahead cross section of equation (8) appears in the ion transport of equation (26), 
only the forward propagating components appear in the transport of neutrons given by 


^jkjor 


{E,E') 


f O' ■/.(E,E',fl,fl')dfl' ^ j^n 

J 47T 


(28) 


In equation (26), the use of the straight-ahead approximation allows Q to be viewed as a parameter. The equation 
may subsequently be solved according to the boundary condition set as the incident fluence from the direction = 
Qo as shown in Fig. 9. The coupling to the isotropically produced neutron field for each Qo is given by 


[J2.V + a^{E)]cl^„^^Jx,^l,E) = dE^ 


(29) 


Once the forward fluence is evaluated for direction Qo over the domain ofx and E, one must yet solve equation (29) 
for the induced isotropic neutron fluence within the shield configuration. Note that the last term in equation (29) is 
an isotropic source of neutrons from collisions of the forward propagating fluence given as 




(30) 


for which the forward fluence i^/, j^g^(x,flg,E') is a function of the penetration depth z (Fig. 9) along the direction 
Qo. This is denoted by introducing Xn E] = where z(x) is the penetration depth 

along Qo to the point x. Equation (29) is rewritten as 


[^l.V + a^{E)]d^„,JxM,E) = £^a^jE,E\Q,Q’)<f>^^^Jx,fl’,E’)dQ’dE’ ^ 

+X„,^J<x),^,^0’E] 

Note that equation (31) also provides a basis for introducing anisotropic source terms as well. Whereas the 
forward propagating solution has its source fixed by the boundary condition, the isotropic solution (having no 
inbound fluence) is driven by internal sources generated by the collisions of the inbound fluence as modified by the 
forward propagator in penetrating to x along direction Qo. In evaluating the penetration depth, one must account for 
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differences in the material composition that determines the attenuation in reaching x from the boundary in addition 
to the strength of the neutron source term at location x. 

The direction Q may be viewed as a parameter in equation (31), allowing the equation to be solved along an 
arbitrary number of directions used to construct a numerical representation of 4>^ (*, fJ, E) . The Q dependence of 

comes from the distance to the boundary t(Q) and the penetration depth z along Qo. The isotropic 
nature of the neutron source of equation (31) effectively decouples the incident direction Qo from the neutron 
direction of propagation Q. To solve equation (31), one must know the neutron fluence, (j)^-^g{x,n,E) , at every x 

from all directions Q' approaching x from adjacent locations, as was the case of equation (14). Although the use of 
the approximation in (15), based on the assumption that losses along the ray Q are compensated by diffusion from 
adjacent rays, is a demonstrated useful assumption in a semi-infinite slab (see Fig. 8 and Wilson et al. [2014], Slaba 
et al. [2010a, 2011b]), it is generally an overestimate in a sphere (or any finite object). Using equation (15) in 
evaluation of equation (31) seems a reasonable approximation in most human rated vehicles, leading to an extension 
of transport along an arbitrary ray Q as a solution to the bi-directional equations 

[Q. V + aJE)]d^^^^Jx,fl,E) = J^a^^j{E,E')d^^^^JxME')dE' 

+ J^<^nn,tiE,E')^^^^Jx-n,E<)dE' (32) 

where the backward moving field satisfies a similar equation obtained by replacing Q by -Q in equation (32) as 

[-Q. V + = J^a^^j{E,E')cl^^^^Jx-n,E')dE< 

+ J^<^nnA^>E')^^^,Jx,n,E')dE' (33) 

These isotropically produced neutrons produce isotropic light ions in collisions with the media resulting in a light 
ion fluence approximately given by 


n»Vcj)^{x,n,E) 


A^dE 


S,{E) 


/ •oo 

^ ^jn(E,E')cP„^^Jx,n.E')dE\ 


(34) 


where (j)^-^g{x,fl,E’) is the isotropically produced neutron fluence. Equations (32) and (33) are solved 

simultaneously for no additional external fluence incident on the boundary of the media directed along Qo (see Fig. 
9). Note that the implicit assumption of equations (32) and (33) is that lateral diffusion from the ray Q is fully 
compensated by diffusion from adjacent rays. This was demonstrated for slab geometry as a reasonable result as 
shown in Figs. 6-8. The lateral leakage for incident particles along Qo in this case is found to a first order 
approximation by solving equations (32) and (33) along each Q, between the ray limits {-t(-Q,), t(Q/)} for a given 
set of A directions Q,. A set of directions from A= 1 to N= 22 was introduced by Wilson et al. [2014], and given in 
Appendix B. Convergence tests are applied in the next section to verify the ability of discrete ray numbers in 
describing the isotropic radiation fields. 

A Simple Inhomogeneous Geometry 

Consider now a 15 g/cm^ radius tissue equivalent sphere (ICRU) shielded by a 20 g/cm^ spherical shell of 
aluminum, as shown in Fig. 10. The penetration depth now generally includes aluminum and tissue thicknesses. The 
solution along any direction Q, must generally include propagation through both media from the outer limits of the 
object along the ray ±Q, to the evaluation point. In this example, the solution within the ICRU sphere serves as 
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proxy for the astronaut (ICRU phantom, see ICRU [1993]). In particular, the D*{td, Qo) and//*(trf, Qo) as defined by 
ICRU are used. With verification of the methods in this simplified geometry, these methods will be ready to apply to 
more realistic astronaut and spacecraft geometry in subsequent studies. Simple geometry is chosen to allow efficient 
MC simulation for comparison with the calculation of quantities using the above formalism. In this way, the solution 
methodology is verified and experience gained in the formalism in preparation for connecting to the general 
geometric descriptions generated by the engineering design process for which MC simulation is difficult to 
implement. Target points (detectors) are placed at depths along the z-axis {x = 0,y = 0). Most results will be shown 
at or near detectors located at the top of the tissue sphere (0 g/cm^) and bottom of the tissue sphere (30 g/cm^) as 
shown in the right pane of Fig. 10. These depths are being shown to highlight transition effects that occur near 
material interfaces, as previously shown in Figs. 6 and 7. 

The external radiation environment (external source boundary condition) is assumed to be anti-parallel with 
the z-axis, uniform in the x-y plane, positioned above the inhomogeneous sphere, and directed down onto the top of 
the sphere as shown in Fig. 10. The added complexity in the calculation is the separate evaluation of the tissue 
thickness and shield thickness along each Q, along with the evaluation of the penetration through the shield and 
tissue along Qo for each Q for use in the transport procedure. Note that the penetration depths for sources in the 
lower hemisphere will generally have double penetration through the shield, once at the top and again at the bottom 
after passing through the tissue equivalent sphere. 


rio 






Fig. 10. Geometric relations for a sphere of tissue shielded by a spherical shell of aluminum. The left pane shows the 
geometry with variables related to ray-tracing description. The right pane shows dimensions of geometry (in g/cm^) 
and locations of detectors where results will be shown and compared. Both geometries drawn in cm scale. 


Field values are to be evaluated at a point x where the thickness of material t is known over an array of 
directions, Q, (i.e. t, = Note that each Q, for i odd is backed by Q,+i = -Q,. To complete the solution 

methodology, the directions of propagation, Q„ must be selected, and the penetration depth, Zp, for a given Qo must 
be evaluated. The N directions, Q,, are chosen such that each ray accounts for the same solid angle, AQ = 4tz/N, on 
the unit sphere. The evaluation will be made for a single direction of incidence, Qo, but this could include an array 
of incident directions for which the solutions would be integrated after all transport procedures have been 
completed. In that low-energy neutron transport is dominated by diffusive processes, it is anticipated that the 
solution will rapidly converge for low values ofN. 

The source of neutrons induced within the inhomogeneous sphere by the Webber SPE event is shown in 
Fig. 11. Note that the excess attenuation in tissue relative to aluminum for the primary protons is seen as the slight 
shadow cast by the tissue sphere, and the weaker neutron source strength within the tissue shifts the neutron 
production to lesser values easily seen at the top of the tissue sphere. The source of neutrons induced by four 
components of the 1977 solar minimum GCR ions of hydrogen, helium, carbon, and iron is shown in Fig. 12. The 
greater penetrating power of the GCR is apparent in the figures resulting from the much greater energy of the GCR 
ions. Also, note that the fluence of the forward propagating light ions increases with penetration depth, forming the 
most intense neutron production near the exit. While the GCR hydrogen and helium result in a shadow behind the 
tissue sphere, the more complex ions (carbon and iron in the present example) are more effectively fragmented on 
the hydrogen contained in the tissue, forming a more intense source of neutrons behind the tissue sphere. 
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Fig. 11. Integral {E > 1 MeV) isotropic neutron source (particles/(g-event)) induced by the Webber SPE incident on 

an ICRU sphere and aluminum spherical shell. 
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Fig. 12. Integral (E > 1 MeV) isotropic neutron source (particles/(g-day)) induced by 1977 solar minimum ions of 


hydrogen, helium, carbon, and iron incident on an ICRU sphere and aluminum spherical shell. 


The resulting neutron fluence depends not only on the distance, t(Q), to the boundary but the intensity of 
the source along the ray in reaching the boundary. Clearly for solar particle events, the most intense regions must be 
adequately represented in the domain of the solution. The convergence rate for increasing number of rays N will be 
tested in the following sections. The number of rays, Q,, for the current example is now discussed. 
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N = 1, Straight-ahead approximation 


For each Qo, only one value of Q, is allowed and is taken to be that of the incident direction Qo. The transport 
distance x along Qi and penetration depth z are the same, reducing complication in the formulation ioxN = 1. In this 
case, leakage from the media is only through the distal face. The structure of this code is quite different than the 
original HZETRN code, although computational results are identical, and the original HZETRN code provides a 
verification test-case of the new 3DHZETRN code. For N = \, the comparisons in the prior sections describing the 
verification and validation of HZETRN are fully applicable. If one restricts interest to I, then the HZETRN code 
is preferred in that computational efficiency is higher, but the added overhead of 3DHZETRN for 1 is minimal. 
In the following sections, comparisons between 3DHZETRN and MC simulations will also include results from the 
HZETRN code (A/^=l) to highlight shortcomings of the straight-ahead approximation and improvements of the 3D 
corrections presented herein. 

N = 2, Bi-directional approximation 

The first step towards a full 3D code is the treatment of backward propagating particles. This was first 
accomplished by Clowdsley et al. [2000, 2001] and later improved by Slaba et al. [2010a] with demonstrated 
improved results [Slaba et al. 2010a, 2011b; Heinbockel et al. 2011a, 2011b]. In this case, propagation of the 
isotropic neutrons is represented in two directions by Qi = Qo and Q 2 = -Qo. Therefore, the penetration depth z is 
equal to x along Qi and equal to t(Qi) + t(Q 2 ) - x along Q 2 . Note that lateral leakage is not represented in this 
approximation, although leakage at the forward and backward boundaries is represented. The neutron transport in 
semi-infinite slab geometry has extensive MC results for verification for which this bi-directional approximation, 
ignoring lateral leakage, is most appropriate as demonstrated in Fig. 8. This approximation is appropriate for slab 
geometry and works well in applications in which the shield geometry is a collection of thin plates. The next step for 
treatment of lateral leakage is now discussed. 

N > 2, N-stream approximation 

For general N, the transport formalism requires the penetration depth, Zp, to be known as a function of 
transport depth, xt, along the propagation direction Q for each Qo. In order to simplify algebraic manipulation and 
expressions, the center of the sphere is assumed to be positioned at the origin with the positive z-axis anti-parallel to 
the source direction Qo, as shown in the left pane of Fig. 10. Due to the simplified geometry and source orientation, 
an analytic expression for the penetration depth through the aluminum shield and tissue target may be determined. 
The point xo may be written as 


X, 


'0 


X — 


Q [ (fj + it) ~ (^r ^ ^ ) 


(35) 


where ts and U are the known shield and target thicknesses along direction -Q from x to the outer boundary, 
respectively. The symbols ^ and ^ represent the transport depth through the shield and target along Q, 

respectively. The points and may be written in terms of Xo as 



(36) 


(37) 


where and are the penetration thicknesses through the aluminum shield and tissue target, respectively. The 

expression for z^^^ is obtained first by writing equation (37) in component form and substituting into the equation of 
a tissue sphere with radius 15 g/cm^. This yields 
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(38) 


= a:„*r2n 


15 ^ 




■'0 ““0 


The expression for is obtained by writing equation (36) in component form and substituting into the equation of 
an aluminum shell with thickness 20 g/cm^ surrounding a tissue sphere with radius 15 g/cm^. This yields 



(39) 


For preliminary testing and benchmark comparisons against MC codes, the solutions for N = 6, 10, 14, 18, 
and 22 are considered separately, in addition to the simple N~ I and N ^2 cases discussed above. The distributions 
for all N are shown in Fig. 13. Note that the N = 6, 14, and 22 cases include rays perpendicular to the z-axis (i.e. 
tangent to the sphere surface at the top and bottom detector locations). On the surface of the sphere, these rays 
correspond to shielding by aluminum only with no intervening tissue. This is corrected in the 10 and 18 
distributions. The directional elements for each N are given in the Appendix B, along with a brief description of how 
the elements were determined. For each N, the omni-directional flux is obtained by integrating the fluxes along each 
ray direction after the transport procedure. 
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Fig. 13. Directional components for 1,2, 6, 



Preliminary testing 

One very large solar particle event occurred on February 23, 1956 and has served as an important event to 
not only space operations but for high altitude aircraft as well [Foelsche et al. 1974]. As a result, this event has 
become a test environment used over many decades and has been (at times) approximated by the Webber spectrum 
given by equation (24). Of interest has been the dose and dose equivalent within the tissue sphere for this Webber 
spectrum. Dosimetric quantities are evaluated at depths in the tissue sphere at td = 0, 0.007, 0.3, 5, 25, 29, 29.7, 
29.993, 30 g/cm^ with results in Tables 3 and 4 for A = 1 through N = 22. The dose even on the surface of the tissue 
sphere is seen to converge rapidly with increasing N. The local fields appear dominated by transport through the 
aluminum shell giving improved convergence behavior over the case concerning the surface of the aluminum shield 
[Wilson et al. 2014]. The rate of convergence of dose equivalent is somewhat slower as the details on the penetration 
of charged particles produced by neutrons are more dependent on angular factors. 

Dose and dose equivalent convergence rates are reasonably rapid for all locations within the tissue sphere 
except near the tissue/aluminum interface as shown in Tables 3 and 4. At each depth, results appear to 
asymptotically approach a common value, as will be further demonstrated for transported fluence spectra. The 
induced particle fluence is shown in Fig. 14 at the near and distal interfaces. Also shown are the lowest order N = \ 
and N ~ 2 solutions of the HZETRN code. The N = I solution has no contributions from the backward leaking 
neutrons that have been moderated within the tissue sphere, as evidenced by the incorrect spectral shape at the 
lowest energies. The incorrect spectral shape is largely corrected by the N = 2 bi-directional approximation as seen 
in Fig. 14. 
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Table 3. Webber SPE dose as a function of number of rays, N, at various depths including those defined by ICRU. 


d (g/cm^) 



Dose (cGy/event) 



1 

N^2 

N=6 

o 

II 

14 

oo 

II 

N=22 

0 

7.01 

7.01 

7.03 

7.03 

7.02 

7.03 

7.03 

0.007 

7.00 

7.00 

7.08 

7.03 

7.05 

7.03 

7.04 

0.3 

6.64 

6.70 

6.73 

6.71 

6.72 

6.72 

6.72 

1 

6.07 

6.12 

6.14 

6.13 

6.13 

6.13 

6.13 

5 

3.83 

3.87 

3.86 

3.86 

3.86 

3.86 

3.86 

25 

0.72 

0.73 

0.71 

0.71 

0.71 

0.71 

0.71 

29 

0.56 

0.56 

0.54 

0.55 

0.55 

0.55 

0.55 

29.7 

0.53 

0.54 

0.52 

0.52 

0.52 

0.52 

0.52 

29.993 

0.52 

0.53 

0.51 

0.51 

0.51 

0.51 

0.51 

30 

0.52 

0.52 

0.50 

0.51 

0.51 

0.51 

0.51 


Table 4. Webber SPE dose equivalent as a function of number of rays, N, at various depths including those defined 
by I CRU. 


d (g/cm^) 



Dose equivalent (cSv/event) 



N= 1 

N=2 

N=6 

N= 10 

N= 14 

oo 

II 

N=22 

0 

11.33 

11.15 

11.61 

11.35 

11.41 

11.45 

11.42 

0.007 

11.18 

11.01 

12.10 

11.31 

11.58 

11.42 

11.51 

0.3 

10.05 

10.45 

10.92 

10.62 

10.71 

10.69 

10.70 

1 

9.25 

9.63 

9.82 

9.68 

9.73 

9.73 

9.73 

5 

6.11 

6.36 

6.20 

6.17 

6.20 

6.20 

6.19 

25 

1.46 

1.50 

1.29 

1.35 

1.34 

1.33 

1.34 

29 

1.19 

1.23 

1.06 

1.11 

1.11 

1.10 

1.10 
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1.15 

1.19 

1.05 

1.08 

1.09 

1.08 

1.07 

29.993 

1.13 

1.16 

1.09 

1.06 

1.09 

1.06 

1.07 

30 

1.13 

1.12 

0.98 

1.01 

1.01 

1.01 

1.00 



Fig. 14. Webber SPE convergence test in 30 g/cm^ diameter tissue sphere within a 20 g/cm^ spherical shell of 

aluminum at the two interfaces. 


The solution methodology appears to be converged for N > 10. This is supported in Tables 3 and 4 by 
comparing exposure values for different N ai a fixed depth. The variation is found to be less than 1% for dose and 
less than 2% for dose equivalent for N> 10 at all depths. Convergence for A^> 10 is also supported in Fig. 14 in that 
spectra for neutrons, protons, and alphas are nearly indistinguishable. The fluence spectra in the neighborhood of the 
near and distal interface are shown in Fig. 15. A rapid transition occurs in passing from aluminum to tissue in the 
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near interface and in passing from tissue to aluminum in the distal interface. The alpha fluence spectra are shown in 
the neighborhood of the near and distal interfaces on a linear scale in Fig. 16. Qualitatively, the spectral 
characteristics are similar to those near the slab interfaces in Fig. 7. For example, the alphas produced in the 
aluminum shield rapidly form a new equilibrium spectrum within the first few millimeters of tissue as observed in 
the near interface. Similarly, at the distal interface, the leakage of alphas from the tissue into the aluminum shield 
quickly forms a new equilibrium spectrum within the first few millimeters of aluminum. 



Fig. 15. Webber SPE induced fluence within the ICRU sphere/aluminum spherical shell neighboring the near (left 
panel) and distal (right panel) tissue/shield interface with 22. 




Fig. 16. Webber induced alpha fluence within the ICRU sphere/aluminum spherical shell neighboring the near (left 
panel) and distal (right panel) tissue/shield interface with 22. 


MC benchmarks will now be used to verify as much of the solution as possible. It is found that the details 
of the secondary alpha spectra for the Webber SPE boundary condition (for example, see Fig. 16) are difficult to 
verify with MC codes in the spherical geometry without excessively large run-times. For the same reason, defining 
the transition regions near the interfaces for neutrons and protons is also difficult. Consequently, the larger features 
are used to verify the solution of the neutron and proton fiuences at the interfaces and results deeper within the tissue 
sphere will be the focus. If the trends in the solution over larger distances are verified, it is expected that the details 
revealed by 3DHZETRN are reasonably correct. The N = 22 will be used in the MC benchmarks discussed in the 
next section. 
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Monte Carlo Benchmarks 


In 2005, a systematic benchmark exercise was administered as a means of verifying various computational 
procedures resulting from the first phase of code development (1996-2004). One component of transport code 
verification uses established benchmarks using standard environments: the Webber SPE spectrum given in equation 
(24) and the 1977 solar minimum GCR environment. The standard configuration of 30 cm semi-infinite slab of 
water shielded in front by 20 g/cm^ of either aluminum or iron utilizing only dose and dose equivalent (ICRP 26 
quality factor) has been used in the past [Wilson et al. 1991]. It was generally concluded that HZETRN agreed with 
the MC codes as well as the various MC codes agreed among themselves [Wilson et al. 2006; Heinbockel et al. 
2011a, b; Slaba et al. 2010a]. There are currently no compelling reasons to attempt to use MC methods in design and 
operational applications. However, as will be further demonstrated below, the MC codes provide an indispensable 
guide in the development of the 3DHEZTRN code, which does provide a useful basis for space radiation protection. 
Attention is now directed towards benchmarking 3DHZETRN in which the added complications of 
finite/inhomogeneous geometry are addressed. 

Webber SPE Benchmark 

The benchmark SPE spectrum is the Webber [1966] approximation of the February 23, 1956 SPE (not an 
accurate representation but a historically useful test spectrum, see [Scott and Alsmiller 1967]). For the present study, 
the induced neutron and proton fiuence spectra were evaluated in the 70 g/cm^ diameter aluminum/tissue sphere by 
three MC codes (Geant4, PHITS, and FLUKA) and 3DHZETRN {N = 22) with the results at 0 and 30 g/cm^ depths 
in the tissue sphere shown in Fig. 17. Monte Carlo results for protons below 1 MeV have been removed since the 
statistical fluctuations were exceedingly large. 

The proton fiuence predictions among the three codes are in reasonable agreement above 1 MeV with very 
large statistical uncertainty below 1 MeV (not shown). The larger discrepancy lies with the neutron spectra where 
disagreement among the MC is on the same order as for 3DHZETRN. The FLUKA and 3DHZETRN results exhibit 
a similar inflection in the neutron spectrum near 10 to 60 MeV that may be related to the use of an early FLUKA 
approximation to the secondary neutron spectrum used in HZETRN [Ranft 1980, Wilson et al. 1988] and carried 
over into 3DHZETRN. Some of the differences in proton spectra on the back surface may be associated with the 
simplified assumptions of the straight-ahead approximation in light ion transport. 




Fig. 17. Particle spectra at detector locations 0 g/cm^ (left pane) and 30 g/cm^ from top of tissue sphere (radius 15 
g/cm^) surrounded by aluminum shell (thickness 20 g/cm^) exposed to the Webber SPE spectrum. 
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A principal difference in the three codes is the computational efficiency. The time required by the codes to 
run this benchmark is given in CPU seconds used and listed in the first column of Table 5. The 3DHZETRN (N ~ 
22) time of ~1 minute allows for the development of affordable means of design optimization software on a 
commercial laptop while Geant4, PHITS, and FLUKA codes require a multiprocessing cluster environment of a 
thousand or more processors to do a single design evaluation and would still be impractical for design optimization. 
It should be noted that source biasing, physics biasing, or other optimization strategies could significantly reduce the 
Monte Carlo runtimes shown in Table 5. However, the fully optimized run times would still be orders of magnitude 
larger than the 3DHZETRN run times and still require a multi-processor cluster computing environment. 


Table 5. Total CPU seconds required for benchmark calculations. Total number of histories ran in Monte Carlo 
simulations are shown in parentheses. 



Webber SPE 

GCR hydrogen 

GCR helium 

GCR carbon 

GCR iron 

3DHZETRN (V=22) 

Geant4 

FLUKA 

PHITS 

69 

2x 10H1>^10") 
1x10*^ (4x10“*) 
3x10^(2x10“*) 

952 

1x10® (5x10“*) 
4x10^(1x10“*) 
3x10^ (4x10^) 

955 

2x10® (3x10“*) 
2x10^ (3x10“*) 
8x10^(1x10’) 

953 

2x10® (1x10’) 
6x10’ (1x10’) 
7x10’ (4x10®) 

954 

3x10® (3x10®) 
2x10’ (1x10®) 
1x10® (1x10®) 


1977 Solar Minimum GCR benchmark 

The GCR benchmarks utilized the Badhwar-O’Neill 2010 model [O’Neill 2010] to represent the 1977 solar 
minimum conditions (solar modulation parameter set to 475 MV). The hydrogen, helium, carbon and iron spectra 
were considered separately as external source boundary conditions. The results for the neutron and proton fluence 
spectra for an incident 1977 solar minimum hydrogen spectrum are shown in Fig. 18 at the interface depths of 0 and 
30 g/cm^. 




Fig. 18. Same as Fig. 17 but for incident solar minimum GCR hydrogen spectrum. 


The corresponding MC results for secondary alpha particles had large statistical uncertainty across the energy 
domain and are not shown. The GCR hydrogen results are qualitatively similar to results obtained for the Webber 
SPE, only now 3DHZETRN agrees more closely with PHITS over the high energy range. The 3DHZETRN proton 
fluence below 100 MeV shows better agreement with Geant4. The MC computation time used for transporting GCR 
hydrogen is slightly shorter than that for the Webber SPE transport in the same geometry mainly due to the slower 
convergence of the MC solution at lower energies. The 3DHZETRN (V = 22) run time of ~15 minutes is very 
manageable whereas the MC codes take several orders of magnitude longer. 
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The neutron, proton and alpha spectra generated by GCR helium ions are shown in Fig. 19. A large 
contribution to the alpha fluence spectrum is from the incident heliutu ions transmitted without suffering a nuclear 
interaction. There is generally better agreement with Geant4 results than those generated by GCR hydrogen. PHITS 
and FLUKA have a tendency toward greater production of neutrons and protons below 1 GeV, whereas Geant4 is in 
better agreement with 3DHZETRN. 



Fig. 19. Same as Fig. 17 but for incident solar minimutu GCR heliutu spectrum. 
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Fig. 21. Same as Fig. 17 but for incident solar minimum GCR iron spectrum. 


Similar results are shown for GCR carbon and iron ions in Figs. 20 and 21. The comparisons are more 
mixed in these cases. The carbon ion results in Fig. 20 exhibit closest agreement with FLUKA and Geant4 below 10 
MeV for neutrons and with PHITS and FLUKA above 1 GeV. A similar result is found for proton spectra. This is 
likely due to differences in carbon breakup cross sections. The alpha results for 3DHZETRN agree more closely 
with PHITS after penetrating the top aluminum shield and after passing through the tissue sphere. The 3DHZETRN 
results for incident iron tend to agree more closely with the FLUKA results. It is clear that the breakup of HZE into 
light fragments requires additional study. 

To gain further understanding of the light ion flux problem, the AZ =0, 1,2 flux induced by GCR carbon and 
iron environments is evaluated as shown in Figs. 22 and 23. Note that AZ refers to the charge difference between the 
primary and secondary particles. The Z = 6 flux contains the surviving carbon primaries while Z = 5, 4 are purely 
secondary fragments. It is especially clear that there are large differences in nuclear fragment cross sections between 
the codes, and that this is in part the cause of differences in the light ion flux in Fig. 20. Similar conclusions are 
drawn from the GCR iron ion induced flux shown in Fig. 23, resolving these differences is of great importance to 
shielding from these heavy ions. 



Fig. 22. Particle spectra for AZ =0, 1, 2 at detector locations 0 g/cm^ (left pane) and 30 g/cm^ from top of tissue 
sphere (radius 15 g/cirf) surrounded by aluminum shell (thickness 20 g/cm^) exposed to the 1977 solar minimum 

GCR carbon spectrum. 
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Fig. 23. Same as Fig. 22 but for incident solar minimum GCR iron spectrum. 


Although it is difficult to favor one code over another on the basis of the veracity of the solutions obtained as 
significant differences in the light ion spectra exist, there are even greater differences in computational time required 
to evaluate fluence spectra even in this simple geometry (inhomogeneous spherical aluminum shell containing a 
tissue sphere). Still, MC methods provide an important test point for developing codes capable of meeting 
operational and especially design requirements. In addition, computational speed has proven to be important in 
spaceflight validation in LEO where the time structure of the environment can be used to test various environmental 
components as was done using ISS and the Liulin instrument [Wilson et al. 2007, Slaba et al. 2011a, 2013]. The 
main limitation on such studies remains the uncertainty in the environmental models (especially for the trapped 
environment), uncertainty in nuclear cross sections, and vehicle geometry and composition. 

Conclusions 

The HZETRN and 3DHZETRN radiation transport codes for space design are still incomplete in both the 
physical description and corresponding implementation into computational procedures. Even so, the current 
development appears as state-of-the-art in computational procedures for spacecraft design. The present version of 
3DHZETRN has some inherent limitations (aside from the simplified geometry) that carry over from the available 
bi-directional methods in which the main share of light ions are still treated in the straight-ahead approximation as 
coupled to the forward propagating neutrons. Only the isotropically produced neutrons are first transported in 3D 
without the coupling to the isotropically produced light ion components. Once the isotropically produced neutrons 
are evaluated, the final coupling to the light ions is implemented (without further nuclear reactions) in a full 3D 
manner. Hence, only a fraction of the light ions have a 3D treatment of transport effects of unquantified magnitude. 
Future work will require extension to complex geometries and treatment of the cross-coupling of 3D neutron/light 
ion transport. Even so, the present expanded treatment of 3D effects shows great promise in the path to a fully 3D 
transport algorithm. 

Further limitations result from the current HZE nuclear database. The current database used in 3DHZETRN 
only accounts for the breakup of the projectile ions without regard to the struck nucleus breakup. Furthermore, the 
projectile fragments are assumed to have the same velocity as the initiating projectile. Adding these neglected target 
fragments and correcting the projectile fragment spectra will result in a more reliable simulation. 
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Appendix A: Monte Carlo simnlations 

In this section, a brief description is given of the MC codes Geant4, FLUKA, and PHITS used to simulate 
particle spectra shown in Figs. 8 and 16 - 23. A discussion of the source, geometry, and simulation setup used for the 
codes is also given. 

Geant4 version 9.4.6 [Agostinelli et al. 2003] was used in the present calculations with the physics list 
QGSP BERT HP. This physics list utilizes the quark gluon string model (QGS) for high energy interactions with 
nucleons, pions, and nuclei. Post-interaction nuclear de-excitation is handled by the precompound (P) model. The 
Bertini cascade (BERT) model is used for interactions below 10 GeV. A high-precision (HP) tracking model is used 
for neutrons below 20 MeV. Nucleus-nucleus collisions are represented by the native Quantum Molecular Dynamics 
(QMD) model [Koi 2008]. More information about the chosen physics list can be found at the Geant4 website 
[Geant4 2012a]. The QGSP BERT HP physics list has been recommended by the Geant4 collaboration for high 
energy physics applications [Geant4 2012b] and has been used, with some slight variation, in other space radiation 
related studies [Bernabeu and Casanova 2007, Hayatsu et al. 2008, Martinez and Kingston 2012, Slaba et al. 2013]. 

FLUKA is a MC program that contains fully integrated physics and performs the transport of elementary 
particles and ions through materials. It has the ability to transport and calculate interactions with all elementary 
hadrons, light and heavy ions, and electrons and photons over an energy range which extends up to 10"^ TeV for all 
particles, and down to thermal energies for neutrons [Battistoni et al. 2007]. The nuclear models are integrated into 
the software for all particles and energies, with the exception of neutrons below 20 MeV, where standard 
international evaluated data files are used. For nucleus-nucleus interactions above 5 GeV/n, FLUKA utilizes 
DPMJET-III with a modified initialization procedure, and for energies between 0. 1 and 5 GeV/n, a modified RQMD 
model is used [Roesler et al. 2001, Sorge et al. 1989a, Sorge et al. 1989b]. For low energy nucleus interactions 
below 0.1 GeV/n, the Boltzmann Master Equation is used [Cerutti et al. 2006]. All models are followed by a 
generalized intra-nuclear cascade process followed by a pre-equilibrium stage and then an evaporation- 
fragmentation-fission stage. This code has been extensively benchmarked against available accelerator and cosmic 
ray experimental data, at beam energies as low as a few MeV and as large as cosmic ray energies [Aarnio et al. 
1993; Andersen et al. 2004]. 

PHITS has been used extensively in space radiation transport. It simulates ionization through transport 
processes and uses the mean free path to determine when collisions occur [Niita et al. 2006, Sihver et al. 2007]. 
Ionization transport can include angle and energy straggling, showing good agreement for Bragg peaks and 
fragmentation tails [Niita et al. 2006]. Tabulated nuclear cross section data are used for neutron interactions below 
20 MeV, and evaluated data are also used for photon transport below 100 GeV and electron/positron transport below 
100 MeV [Sato et al. 2013]. Medium energy neutron and light ion induced nuclear interactions are simulated 
through the Intra-Nuclear Cascade model of Liege [Boudard et al. 2013], while high energy (>3 GeV) nuclear 
interactions involving neutrons and other hadrons are simulated through the Jet AA Microscopic transport model 
(JAM) [Sato et al. 2013]. The Japan Atomic Energy Research Institute (JAERI) Quantum Molecular Dynamics 
model (JQMD) addresses nuclear interactions induced by deuterons, tritons, helions, and alphas at energies above 3 
GeV/n and heavier charged particle induced nuclear interactions from 10 MeV/n to 100 GeV/n [Sato et al. 2013]. 
Once the nuclear interaction dynamic portion is complete, the General Evaporation Model is used to address nuclear 
de-excitation through particle emission in a statistical manner [Niita et al. 2006, Sihver et al. 2007, Sato et al. 2013]. 
The default total nucleus-nucleus cross sections used in PHITS are the same as those that have been used in versions 
of HZETRN [Niita et al. 2006, Sato et al. 2013]. 

In Fig. 8, the geometry setup was for a 30 g/cm^ ICRU tissue slab shielded in front and back by 20 g/cm^ of 
aluminum. All slabs are assumed to have infinite lateral dimensions. Particle spectra were tallied through surfaces at 
the depths indicated in the figure. In Figs. 17 - 21, the geometry setup was for an ICRU tissue sphere with radius 15 
g/cm^ surrounded by an aluminum shell with thickness 20 g/cm^. The external source was applied uniformly onto 
the top of the sphere (covering the entire sphere) parallel with the z-axis as shown in Fig. 10. Particle spectra were 
tallied at spherical tissue detectors (radius 2.5 mm) placed along the z-axis down through the sphere. Convergence 
tests on the detector radius were performed to arrive at the final choice. In general, the chosen detector radius is 
large enough to allow reasonable tally statistics, but small enough to effectively prevent fragmentation events from 
occurring within the detector. 

For all of the GCR simulations, the source energy spectrum was sampled directly from the input spectrum 
(i.e. no source energy biasing). For the SPE simulations, the source energy spectrum was biased so that source 
sampling was performed from a uniform distribution covering 0.01 MeV to 2500 MeV. This biasing technique 
places greater emphasis on the higher energy portion of the SPE spectrum improving tally statistics for secondary 
neutrons. 
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In Figs. 22 and 23, the simulation setup is identical to the other GCR simulations, except for the following 
changes to improve efficiency and statistical uncertainty. First, all particles produced in the simulation with AZ > 2 
(i.e. anything lighter than Be for the GCR carbon benchmark and anything lighter than Cr for GCR iron benchmark) 
were terminated. Second, instead of covering the entire geometry, the source irradiation plane was restricted to a 
radius of 5 cm. Both modifications greatly reduced simulation run time allowing more histories to be simulated, 
thereby reducing statistical uncertainties. The first modification is justified since HZE particles produced within 
shielding receive little or no contributions from lighter mass particles [Wilson and Badavi 1986, Wilson et al. 1991]. 
The exception to this occurs at low energies, less than a few MeV/n, which receive contributions from nucleon- 
induced target fragmentation. These energies are not shown in Figs. 22 and 23, and MC codes generally have 
difficulty resolving this part of the spectrum due to the short residual ranges. The second modification is justified 
since HZE particles are highly peaked in the forward direction [Wilson and Badavi 1986, Wilson et al. 1991]. 
Convergence testing with Geant4 revealed that no particles were scored in a detector (out of 1 billion primaries) if 
the primary ion struck the geometry with a radius greater than 2 cm. The larger source radius of 5 cm was selected to 
ensure even the largest production angles in HZE collisions were accounted. 

Appendix B: Evaluation of Q, for W = 6, 10, 14, 18, 22 

The Q; are herein evaluated with Qo as the top hemispheric pole so that Qi = Qo, while the bottom pole is at 
02 = -Oo and are the central rays of the two polar regions. The remaining directions denoted by Q/ for i = 3...N are 
symmetrically distributed about the mid-latitudes. Generally in application, the 0,’s must be rotated to the body 
coordinates of the transport media. This is accomplished by a rotation matrix R(fti, ^o) determined from the polar 
vector Oo relative to the body coordinate system (note, 6 denotes co-latitude with values on 0 to n, and (j) denotes 
longitude running from 0 to 2%). The regional boundaries are chosen such that 

r dO, = 47 t / A , (40) 

JAn. ‘ 

where AO, denotes the f*' region, N is the number of regions and is related to the number of Ne segments of Acos^, 
and the number of segments A(j)i plus the two polar regions such that 

N=N^N,+2, (41) 

and defined such that 

N N 

X]AO, = (42) 

2 = 1 2 = 1 

The polar caps are bound by 


cos9^^l-2/N, (43) 

COS02 = -1 + 2 / A , (44) 

and the corresponding central rays Oi and O 2 lie along the ±z axis. The middle latitudes are divided into N(jNe 
regions extending from the polar cap to the equator in each hemisphere. The corresponding longitudinal boundaries 
are set by the such that the 


A0, =27t/A^. 

Note that is generally greater than or equal to 4. The remaining boundaries in latitude are given as 


( 45 ) 
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(46) 


A cos 9. — (cos 9^ — cos 9^) / Ng. 

For N = 6, there is only one midpoint in latitude taken as 

cos 9^ — (cos 9-^ + cos 9^) / 2 — Q , (47) 

for i>2. The four longitudinal points are given as 

0 = {O,^/2,7r,37r/2}. (48) 

The directional elements for each solid angle region are given by 

rj, = |cos 4>^ sin 9^ , sin 0. sin 9^ , cos 0; (49) 

and are given for N ^ 6 in Table B. 1 . Similar results for other N used in this paper are given in Tables B.2 - B.5. 

Table B.l Directional elements of the Q/ for N~ 6. 


i 

1 

2 

3 

4 

5 

6 

Cx 

0 

0 

1 

-1 

0 

0 

Cy 

0 

0 

0 

0 

1 

-1 

Cz 

1 

-1 

0 

0 

0 

0 


Table B.2 Directional elements of the fl/ for 10. 

i 1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

6x 0 

0 

0.9165 

-0.9165 

0 

0 

-0.9165 

0.9165 

0 

0 

Cy 0 

0 

0 

0 

0.9165 

-0.9165 

0 

0 

-0.9165 

0.9165 

ez 1 

-1 

0.4 

-0.4 

0.4 

-0.4 

0.4 

-0.4 

0.4 

-0.4 


Table B.3 Directional elements of the Q/ for 14. 


i 

1 

3 

5 

7 

9 

11 

13 

Cx 

0 

0.8207 

0 

-0.8207 

0 

1 

0 

Cy 

0 

0 

0.8207 

0 

-0.8207 

0 

1 

Cz 

1 

0.5714 

0.5714 

0.5714 

0.5714 

0 

0 

i 

2 

4 

6 

8 

10 

12 

14 

Cx 

0 

-0.8207 

0 

0.8207 

0 

-1 

0 

Cy 

0 

0 

-0.8207 

0 

0.8207 

0 

-1 

Cz 

-1 

-0.5714 

-0.5714 

-0.5714 

-0.5714 

0 

0 


Table B.4 Directional elements of the Q/ for 18. 


i 

1 

3 

5 

7 

9 

11 

13 

15 

17 

Cx 

0 

0.7454 

0 

-0.7454 

0 

0.9750 

0 

-0.9750 

0 

Cy 

0 

0 

0.7454 

0 

-0.7454 

0 

0.9750 

0 

-0.9750 

Cz 

1 

0.6667 

0.6667 

0.6667 

0.6667 

0.2222 

0.2222 

0.2222 

0.2222 

i 

2 

4 

6 

8 

10 

12 

14 

16 

18 

Cx 

0 

-0.7454 

0 

0.7454 

0 

-0.9750 

0 

0.9750 

0 

Cy 

0 

0 

-0.7454 

0 

0.7454 

0 

-0.9750 

0 

0.9750 

Cz 

-1 

-0.6667 

-0.6667 

-0.6667 

-0.6667 

-0.2222 

-0.2222 

-0.2222 

-0.2222 
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Table B.5 Directional elements of the Q/ for 22. 


i 

1 

3 

5 

7 

9 

11 

13 

15 

17 

19 

21 

ex 

0 

0.6863 

0 

-0.6863 

0 

0.9315 

0 

-0.9315 

0 

1 

0 

ey 

0 

0 

0.6863 

0 

-0.6863 

0 

0.9315 

0 

-0.9315 

0 

1 

Cz 

1 

0.7273 

0.7273 

0.7273 

0.7273 

0.3636 

0.3636 

0.3636 

0.3636 

0 

0 

i 

2 

4 

6 

8 

10 

12 

14 

16 

18 

20 

22 

Cx 

0 

-0.6863 

0 

0.6863 

0 

-0.9315 

0 

0.9315 

0 

-1 

0 

ey 

0 

0 

-0.6863 

0 

0.6863 

0 

-0.9315 

0 

0.9315 

0 

-1 

ez 

-1 

-0.7273 

-0.7273 

-0.7273 

-0.7273 

-0.3636 

-0.3636 

-0.3636 

-0.3636 

0 

0 
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